#include "fermiqcd.h"

# smc@physics.syr.edu

void set_repr(gauge_field &V, gauge_field &U) {
  mdp_site x(U.lattice());
  int nc=U.nc;
  const int nd=10;
  static vector<mdp_matrix> S(nd);
  static bool done=false;
  if(!done) {
	  done=true;
	  for(int i=0; i<nd; i++) S[i]=mdp_identity(nc); // edit this line!
  }
  forallsites(x)
  	for(int a=0; a<nd; a++)
  		for(int b=0; b<nd; b++)
  			for(int mu=0; mu<U.ndim; mu++)
	  			V(x,mu,a,b)=trace(S[a]*U(x,mu)*S[b]*transpose(U(x,mu)));
}

void set_repr_optimized(gauge_field &V, gauge_field &U) {
  mdp_site x(U.lattice());
  int nc=U.nc;
  const int nd=10;
  vector<mdp_matrix> P(nd),Q(nd);
  static vector<mdp_matrix> S(nd);
  static bool done=false;
  if(!done) {
	  done=true;
	  for(int i=0; i<nd; i++) {
		  S[i]=mdp_identity(nc); // edit this line!
		  P[i]=mdp_zero(nc);
		  Q[i]=mdp_zero(nc);
  	   }
  }
  forallsites(x)
    for(int mu=0; mu<U.ndim; mu++) {
  		for(int a=0; a<nd; a++) {
 			for(int i=0; i<nc; i++)
 				for(int j=0; j<nc; j++) {
					P[a](i,j)=0;
					Q[a](i,j)=0;
		 			for(int k=0; k<nc; k++) {
						P[a](i,j)+=S[a](i,k)*U(x,mu,k,j); // S*U
						Q[a](i,j)+=S[a](i,k)*U(x,mu,j,k); // S*U^t
					}
				}
		}
	  	for(int a=0; a<nd; a++)
	  		for(int b=0; b<nd; b++)
	  			for(int i=0; i<nc; i++)
	  				for(int j=0; j<nc; j++)
		  			V(x,mu,a,b)=P[a](i,j)*Q[b](j,i);
	}
}

int main(int argc, char** argv) {
  mdp.open_wormholes(argc,argv);
  define_base_matrices("FERMILAB");
  int nc=3;
  int L[]={10,10,10,10};
  mdp_lattice lattice(4,L);
  gauge_field U(lattice,nc);
  fermi_field psi(lattice,5);
  fermi_field chi(lattice,5);
  gauge_field V(lattice,5);
  coefficients quark;
  quark["kappa"]=.1;
  quark["c_{SW}"]=0.0;

  U.load(argv[1]);
  psi.load(argv[2]);

  set_repr(V,U); ///

  mul_invQ(chi,psi,V,quark,1e-8);

  chi.save(argv[3]);

  mdp.close_wormholes();
  return 0;
}
